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AN IRREVERSIBLE CONSTITUTIVE LAW FOR MODELING THE DELAMINATION 

PROCESS USING INTERFACE ELEMENTS 


VINAY K. GOYAL*, ERIC R. JOHNSON*, CARLOS G. DAVILA*, AND NAVIN JAUNKY§ 

Abstract. An irreversible constitutive law is postulated for the formulation of interface elements to 
predict initiation and progression of delamination in composite structures. An exponential function is used 
for the constitutive law such that it satisfies a multi-axial stress criterion for the onset of delamination, and 
satisfies a mixed mode fracture criterion for the progression of delamination. A damage parameter is included 
to prevent the restoration of the previous cohesive state between the interfacial surfaces. To demonstrate 
the irreversibility capability of the constitutive law, steady-state crack growth is simulated for quasi-static 
loading-unloading cycle of various fracture test specimens. 

Key words, composite structures, progressive failure, ply damage mode, intradamage mode, interlami- 
nar damage mode, delamination, interface elements, decohesion elements, buckling, postbuckling 

Subject classification. Structural Mechanics 

1. Introduction. Delamination in composite structures usually originates from geometric discontinu- 
ities and material defects such as free edges, dropped plies, re-entrant corners, notches, and transverse matrix 
cracks. Recently, significant progress has been made in the development of tools to predict intralaminar dam- 
age, which often precedes the onset of delamination. Delamination can be a major failure mode in composites 
structures and can lead to significant loss of structural integrity. Fracture mechanics techniques such as the 
virtual crack closure technique (VCCT)[1, 2] has been successfully used in the prediction of delamination 
growth. However, an initial delaminated area must be predefined and a self-similar crack growth is assumed. 

To overcome the limitations associated with the VCCT, interface elements can be located between 
composite lamina to simulate initiation of delamination and non-self-similar growth of delamination cracks 
without specifying an initial crack. Delamination is initiated when the interlaminar traction attains the 
maximum interfacial strength, and the delamination front is advanced when the local surface fracture energy 
is consumed. A softening constitutive law that relates tractions to the relative displacements is generally 
used to formulate interface elements. The softening constitutive law is based on the procedures used by 
Dugdale[3] and Barenblatt[4] to find the extent of the plastic zone ahead of the crack tip in ductile metals; 
the size of which is chosen such that the stress singularity from linear elastic fracture mechanics disappears. 
The softening portion of the constitutive law accounts for the complex mechanisms occurring in the volume 
of material ahead of the crack tip by which large amounts of energy are absorbed in the fracture process. 
For laminated composites these mechanisms include nucleation, growth and coalescence of microcavities. 
Hilleborg[5] developed the first comprehensive interface finite element model and applied this method in 
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concrete cracking. Later, Needleman[6] developed a cohesive-decohesive formulation to simulate dynamic 
crack growth in isotropic elastic solids. 

The exact mathematical form of the interfacial constitutive law is less important than its capability to 
represent the maximum interfacial strength and critical fracture energy. Hence, the functional dependence of 
the traction field on the displacement jump across the interface is not uniquely defined. However, functions 
with continuous derivatives have a numerical advantage over functions with discontinuous derivatives when 
used with Newton-Raphson method because the tangent stiffness is smooth. A smooth tangent stiffness 
as a function of the relative opening displacement has been found to mitigate the numerical oscillations 
encountered in using a softening constitutive relation and to eliminate oscillatory convergence difficulties [7], 

The exponential function for the softening constitutive law is smooth and mimics the physics involved 
in the separation of two atoms initially bonded [8]. This form of the constitutive law has been used in the 
analysis of crack initiation, dynamic growth, branching, and arrest in homogeneous materials [9] . Shahwan 
and Waas[10] used it to study delamination of composite structures caused by compressive loads. The various 
exponential constitutive laws that have been successfully employed to simulate delamination are based on 
the assumption that the consumed local surface fracture energy can be recovered. This assumption is invalid 
for structural systems with stresses that may internally redistribute upon external loading in such a way 
that cracks arrest and cracks faces close. Ortiz and Pandolfi[ll] postulated a damage model and used an 
exponential constitutive law to account for such irreversibilities. A limitation of this model is that the critical 
energy release rates and the maximum interfacial strengths associated with Mode I, Mode II, and Mode III 
fracture cannot be specified separately. 

The present work aims at the establishment of an exponential softening constitutive law that satisfies 
empirical mixed-mode delamination failure criteria for the onset and progression of delamination. An internal 
state variable is included in the constitutive law to permanently damage the internal surfaces that have 
exceeded maximum strength during the deformation process. The paper is structured as follows: (i) mixed- 
mode fracture criteria, (i) mechanics of idealized interface material, (ii) interface finite element, (iii) finite 
element results, and (iv) concluding remarks. 

2. Mixed- Mode Failure Criteria. A quadratic failure criterion based on interlaminar tractions has 
been used to predict onset of delamination [12, 13], To simulate the progression of delamination under mixed- 
mode loading conditions, the power law form of the fracture criterion that includes Mode I, Mode II and Mode 
III interaction has been successfully used with a bilinear constitutive law[14, 15, 16], Davila and Camanho[17] 
developed a bilinear constitutive law that can be used with any mixed-mode failure criterion. To the authors’ 
knowledge, no work has been found incorporating an empirical failure criteria into the exponential softening 
constitutive law. A brief description of the failure criteria used in this paper is presented next. 

2.1. Criterion for the Onset of Delamination. Under pure Mode I, Mode II, or pure Mode III 
loading, the onset of delamination occurs when the corresponding interlaminar traction exceeds its respective 
maximum interfacial strength. However, under mixed-mode loading, delamination onset may occur before 
each traction component reaches its maximum interfacial strength. An expression that considers the inter- 
action between the traction components under mixed- mode loading is the multi-axial stress criterion given 
as 
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where 7', is the interlaminar traction component associated with the j-direction, T? is the maximum in- 
terlaminar traction, and (T 3 ) = T 3 , if T 3 > 0, otherwise it is zero. This function has been included to 
emphasize that the normal compressive traction T 3 does not contribute to the onset of delamination. In 
Equation 2.1, T e is an effective normalized traction, and a > 2 is a real number that determines the shape 
of the tri-dimensional failure surface. Delamination does not initiate if T e < 1, and initiates when T e = 1. 
With a =* 2, one recovers the quadratic delamination interaction. The failure surface for a = 2 is a convex 
semi-sphere in the space of normalized tractions Tj /7J , j = 1 , 2, 3. As the value of a is increased, the failure 
surface approaches a half-cube surface. 


2.2. Criterion for Progression of Delamination. Delamination propagates when the energy release 
rate equals its critical value under pure Mode I, Mode II, or pure Mode III fracture. Generally, delamination 
growth occurs under mixed-mode loading. Under this type of loading, delamination growth might occur 
before any of the energy release rate components attains its individual critical value. The power law criterion 
based on the one proposed by Whitcomb[18] is 


(2.2) 
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where Gj is the energy release rate under Mode j fracture, and Gj C is the single-mode critical energy release 
rate for j = I, II, III. The material parameter a defines the shape of the failure locus. For a = 2, one 
recovers the linear interaction criterion[19]. The shape of the failure locus is a triangular surface. The 
shape of the failure surface approaches a 1/8-cube surface as a increases from 2. Reeder [20] evaluated 
different fracture criteria for mixed- mode delamination in a brittle graphite/epoxy composite, a toughened 
graphite/epoxy composite, and a tough graphite/ thermoplastic composite using the mixed-mode bending 
(MMB) test specimen. The power law criterion was a reasonable fit to the test data for the three different 
materials. Thus, the failure criterion in Equation 2.2 is incorporated into the constitutive law of the interface 
material. 


3. Mechanics of the Idealized Interface Material. The interfacial material is bounded by an upper 
surface S + and lower surface S~ . These surfaces are assumed to coincide with a reference surface S° in 
the undeformed configuration as is shown in Figure 6.1. Thus, it is said that the interface material is of 
zero thickness. The surfaces S± independently displace and stretch, and are connected by a continuous 
distribution of nonlinear springs that act to resist the Mode I opening or Mode II and Mode III sliding of 
the upper and lower surface. 

It is convenient to define a mid-surface S m where the tractions and relative displacements are evaluated. 
For this purpose, let us consider any two points and P contained in S+ and S~ , which are coincident 
in the undeformed configuration. The locus of the midpoints P m of the line joining P + and P define the 
mid-surface S m of the interface material. Refer to Figure 6.2. The normal and tangential components of 
the traction and relative displacement vector are determined by the local orientation of the mid-surface S m . 
The virtual work done by the cohesive-decohesive tractions is given by 

(3.1) 5W int = J J 5Aj Tj dS m 

for any kinematically admissible relative displacements Aj, where Tj are the interlaminar traction compo- 
nents acting on a unit deformed area conjugate to the relative displacements, and S m is the surface area. 
The resistive tractions that are associated to the relative displacements at the point P m are shown in Figure 
6.2. The interlaminar normal traction is denoted T 3 and the tangential tractions are denoted 7\ and TV 
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In the next section, the components of the relative displacements are obtained in terms of the displace- 
ment field with respect to the undeformed configuration. Next, the constitutive equations that relate the 
relative displacements to the traction field are presented. The kinematics and the constitutive modeling 
describe the mechanics of interface debonding. 

3.1. Kinematics of the Interface Material. The fundamental problem introduced by the interface 
material is the question of how to express the virtual relative displacements between the surfaces in 
terms of virtual displacements. As shown in Figure 6.1, consider a three-dimensional space with Cartesian 
coordinates Xi,i = 1,2,3, and let there be surfaces coincident with S° defined in this space by Xi = 
Xi (? 7 i , 772), where r]i,r]2 are curvilinear coordinates on the surface S° . 

Let the Cartesian coordinates xf 1 = xf (rji, 770), * = 1, 2, 3 describe motion of the upper and lower surfaces 
in the deformed configuration. Any point on in the deformed configuration is related to the same 
point on S° through 

(3.2) xf = Xi + Uf 


where Uf are displacement quantities with respect to the fixed Cartesian coordinate system. The coordinates 
x™ = x™(rji, 172 ), i = 1, 2, 3 define the mid-surface S m given by 

(3-3) X™ = \ (4 + xT) '■ - Xi | i(tV | Ur) 

The surface S m is coincident with S° in the undeformed configuration. As mentioned earlier, the components 
of the relative displacement vector are evaluated at the mid-surface S m . Therefore, the local orientation of 
normal and tangential unit vectors to the surface S m is required. This is, 


(3.4) 


rdxy* <9x£* dxg 1 } T 
l dr]! ’ dr]! ’ dru J 
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where Rji are components of the rotation matrix. Since x” 1 depends on the displacements Uf 1 , the rotation 
matrix also depends on Uf 1 . Therefore, the virtual relative displacement are expressed in terms of the virtual 
displacements as follows, 

(3.11) - 5U7) 

Equation 3.11 is substituted into Equation 3.1 to obtain the expression of the internal virtual work in terms 
of the virtual displacements. This form of the internal virtual work is convenient for the finite element 
formulation. In addition, the differential surface area of the mid-surface dS m in the deformed configuration 
is expressed in the form, 

(3.12) dS m = MdS° 

where dS° is the differential undeformed surface area, and M is a function of the displacement field Ux . 

3.2. Constitutive Equations for the Interface Material. The stress singularities at the crack-tip 
in the linear elasticity solutions, stemming from the sharp slit approximation, cannot be reconciled with any 
realistic local rupture process. From the molecular theory of strength it is known that there exists stress 
limits for which molecular bond rupture occurs. The softening-type of cohesive zone model is intended to 
represent the degradation of the material ahead of the crack-tip. It captures strength-based bond weakening, 
and fracture-based bond rupture. The mechanics of the delamination process comprises three interrelated 
phases: (i) the initiation of delamination, (ii) the evolution of the degradation zone, (iii) and the delamination 
growth. The first phase that takes place is the initiation of delamination, and it is based on a stress limit 
determined experimentally. A stress measure that is used as the limiting value, may involve an interaction of 
interlaminar stresses such as that in Equation 2.1. The second event is the development of a zone ahead of 
the crack-tip that experiences intense deformation such as plastic deformation in metals, elongated voids that 
contains a fibrous structure bridging the crack faces in polymers (crazing), and high density of tiny cracks in 
brittle ceramics. The molecular bonds are weakened and the nonlinear softening behavior is confined in this 
degradation zone, or process zone. The third event, is the growth of delamination, bond- rupture, and it is 
based on a fracture criteria such as Equation 2.2. The constitutive equations to be developed, mathematically 
describe these three delamination phases. The focus of this section is to develop the constitutive equation 
for single-bond rupture based on continuum damage mechanics approach. This particular case is extended 
to mixed-mode delamination. The constitutive equations that are postulated in this section, are shown to 
satisfy the failure criteria for initiation and progression of delamination presented in the previous section. 

Assume that the two material points P+ and P contained in S + and S~ as shown in Figure 6.2 
are connected with a spring. The points are coincident when the spring is unstretched, and a high initial 
spring stiffness restrains the relative displacement of these material points. Under isothermal conditions, the 
traction T that acts to resist the relative displacement of these material points is due to stretching A of the 
spring, and is expressed as 

/ 1 - A? \ 

(3.13) T(A)=T c Aexp ( — - — j 

where A = A/A c , and T c is the maximum bonding strength that occurs at the critical stretching value A c . 
The parameter /? with (3 > 1 and f3 £ R + defines the stretching range for which the bond is weakened before 
complete rupture occurs. It is in this range, that damage accumulates. In Figure 6.3, the traction-stretching 
curve is shown for different values of the parameter (3. The work of debonding per unit area, G c , is given by 
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the area under the traction-stretching curve, or, 


/*00 

(3.14) G c = T(A)dA 

Jo 

= t c A c pP-fl/Pr | exp (^) 

where T[z] is the Euler gamma function of argument z. By prescribing T c , G c , and (3 in Equation 3.14, 
the parameter A c can be computed. The exponential function in Equation 3.13 is a suitable representation 
of a softening constitutive law because with increasing stretching of the spring A, the traction T increases 
to a peak value T c and then decreases until complete debonding occurs. Equation 3.13 is only valid for 
monotonically increasing separation because the work due to stretching is recoverable on unloading. 

An internal state variable d that tracks the damage state of the spring needs to be included in Equation 
3.13 to account for irreversible effects. In the following irreversible law an elastic damage model instead of 
a plastic damage model is assumed, 



Within the framework of continuum damage mechanics, it is possible to impose restrictions on d. It must 
increase as a function of time because thermodynamics requires that the irreversible dissipation associated 
with the debonding process remains semi-positive, i.e., d > 0. An equivalent mathematical expression at a 
discrete time f, is 

(3.16) d (ti) = max (l, d {u ~'\ a£. } ) , d (0) = 1 

with ti > ti- 1 . If the spring is assumed undamaged at t = to, then the initial condition is d^°> = 1. Equation 
3.15 is equivalent to Equation 3.13 if no damage occurs, d = 1, or for monotonic increasing loading, d = A ^ y 
Unloading does not occur linearly to the origin, but with an exponential form. The energy of dissipation 
associated to fatigue is neglected in this work. This assumption is valid in the case of a spring that undergoes 
a small number of loading-unloading cycles. 

Equations 3.15 and 3.16 with (3=1 are used for the traction-stretching curve shown in Figure 6.4. 
The points on the curve labeled 1,..., 6 in this figure, represent a possible loading-unloading cycle of the 
spring connecting terminal points PA The spring is unstretched at point 1. With increasing stretching 
cohesive traction develops to resist the separation. Point 2 is in the stiffening portion of the constitutive 
law, and the action of the spring prevents much separation of points P^ . The onset of delamination occurs 
at point 3, where the traction attains its maximum value. As the spring is stretched beyond the onset of 
delamination to point 4, damage accumulates in the spring and the traction gradually decreases. Assume 
unloading commences at point 4, then the spring begins to contract to point 5. If loading commences from 
point 5, the spring is stretched again to point 6 on the original softening portion of the constitutive law. 
Continued loading from point 6 is along the original softening portion of the constitutive law. The traction 
eventually vanishes as the spring is stretched substantially beyond point 6. 

Equations 3.15 and 3.16 are extended to the mixed- mode delamination case. To develop the constitutive 
equations, it is convenient to normalize the relative displacements A j and the tractions Tj with respect to 
the critical separation values Atj and the maximum interfacial strengths TJ. That is, 

(3.17) Aj = Aj/Aj, % = Tj/T? 
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In reference to Figure 6.2, the components of the normalized relative displacements between with respect 
to the orientation of the surface S m at a point P m is, 

(3.18) v = Ai?i + A 2 r 2 + A 3 r 3 


where fi,r 2,?3 are the unit vectors tangent and normal to the surface S m at a point P m . An effective 
relative displacement A is defined by the norm of v 

(3.19) A — + A| + A| 

We assume that the normalized scalar traction T v acts along the direction of v to resist the effective relative 
displacement A . The proposed constitutive law for the interface material is defined along v, 


(3.20) 


T v { Ai, A 2 , A 3 ) — A Q( Ai, A 2 , A 3 ) 


where Q is a decreasing function of any of the normalized relative displacements A j, j = 1,2,3. The 
components of the traction acting along v, normal and tangent to the mid-surface S m at a point P m is 

(3.21) Tj = T v r j ■ - — - = A j Q 

for j = 1,2,3. The function Q is chosen to satisfy the multi-axial stress criterion in Equation 2.1 for the 
onset of delamination and the mixed-mode fracture criterion in Equation 2.2 and is given by 


(3.22) 


Q = exp 


'2- ^/d-d 

P . 


In the latter equation, the scalar mixed-mode parameter /i is defined such that it couples the normalized 
relative displacements for the opening and sliding modes; i.e., 

(3 23) /i a \ a , 1 a \ a , / A \a\ 1 / a 


( | A 1 T A 2 + (A 3 )“) 


where |Aj| is the absolute value function, and (Ai) = Ai if Ai > 0, otherwise it is zero. The material 
parameter a defines the shape of the failure surface for the onset and progression of delamination. The 
internal state variable d is given by, 

(3.24) = max ^1, Mp)) > = 1 

The constitutive equations are slightly modified to take into consideration the mechanical behavior of 
the interface if penetration of surfaces S ^ is detected. Under these contact conditions the surfaces are 
assumed smooth so that frictional effects can be neglected. When contact is formed between two smooth 
surfaces of adjacent lamina, the equilibrium largely depends upon the distribution of elastic forces in the 
contacting lamina. Adjacent lamina surfaces are under contact at a point P m , P m G S m if the relative 
displacement A 3 between P^ is less than zero. For A 3 < 0, a large repulsive traction T 3 develops to avoid 
interpenetration of the surfaces at P m . From Equations 3.21 to 3.24 we obtain the final form of the 
mixed-mode constitutive equations as 


(3.25) 


f T\ ' 


r a. 

t 2 

> = < 

a 2 

l t 3 J 


, (As) . 


} exp 


2 — — d 

P . 


+ 


0 

0 


> exp 


-<-a 3 > J 


1 + K |A 3 

T~ 
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where k is an interpenetration factor to magnify the repulsive force T 3 , and it is chosen arbitrarily in the 
range k > 1. Equations 3.25 and 3.24 reduce to Equations 3.15 and 3.16, respectively, for single-mode 
delamination. 

The empirical parameters governing the constitutive equations in 3.25 are the critical energy release rates 
G/ c , Giic, Guic, the maximum interfacial strengths Tf, Tf, T^\ the critical separation values Af, A§, A§; 
and the parameter /3. These may be specified based on atomistic models of separation or on a phenomeno- 
logical basis depending whether the separation process is governed by ductile void coalescence or a brittle 
cleavage mechanism. By specifying the critical energy release rates and the maximum interfacial strengths, 
one can obtain the critical separation values. The path independent J-integral along a boundary that con- 
tains the interface material can be used to show that the area under the traction versus separation curve is 
the work of fracture per unit area. Equation 3.14 under pure Mode I, Mode II, or Mode III fracture, is used 
to obtain the critical separation values A^, j = 1,2,3. 


Proof. The exponential constitutive law as specified by Equations 3.24 an d 3.25 satisfy Equation 2.1 for the 
onset of delamination, and Equation 2.2 for the progression of delamination. 


For simplicity, monotonically increasing loading is assumed, i.e., d = max (1,/F 3 ). The effect of inter- 
penetration is also neglected, A 3 > 0. For the onset of delamination, the components of the traction vector 
in Equation 3.25 are substituted into Equation 2.1 to obtain the effective traction T e as 


(3.26) 


T e = ( (A? + A“ + A“) exp 
T -pP 


1 -// 


P 


1/ a. 


= fi exp 


p 


This equation is analogous to Equation 3.13 for single-mode delamination. At the initiation of delamination 
we set T e = 1 in Equation 3.26, to find the only solution which is pP = 1. Hence, the relative displacement 
p = 1 at the onset of delamination when T e = 1. Therefore, the proposed constitutive law, Equation 3.25, 
satisfies Equation 2.1 for the initiation of delamination. 

For the progression of delamination, proportional straining is assumed. The relative displacement asso- 
ciated to the sliding Mode II and Mode III are written as Ai = £2^3 and A 2 = £3^3 with £2 and £3 fixed 
during the loading history. The terms in Equation 2.2 are evaluated as follows, 


G j_ 

G Ic 


a/2 


( f<f* 3 ^MAi, A 2 , A 3 )dA 3 \ 
V f o c °T 3 (0,0,A 3 )dA 3 J 


a/2 


,(l+£ 2 “ + ^) 


2/a 


+ 0l(A 3 ) 


Gn \ 
Giic) 


a/2 


( Jo 1 Ti(Ai, A 2 , A 3 )dA! \ 

V / 0 oo r 1 (A 1 ,0,0)dA 1 J 


e 2 


a/2 


,(l+C 2 “+e 3 “) 


2 A 


+ 02(A 3 ) 
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Ill 


IIIc 


a/2 


( fo " " ^(Ai, A 2 , A 3 )dA 3 \ 
V f o °°T 2 (0,A 2 ,0)dA 2 ) 


tl 


a/2 


,(i + e 2 “ + e 3 “) 


2 /a 


+ 0 3 (A 3 ) 


where <pj( A 3 ), j =e 1, 2, 3 are exponential decaying functions with increasing A 3 . The progression of delam- 
ination occurs when the functions ^(A 3 ) , j = 1,2,3 are virtually zero. Adding the last three equations 
shows that the power criterion in Equation 2.2 is satisfied. 


4. Interface Finite Element. The formulation for the interface element is based on the work of Beer 
[21], An iterative solution procedure is necessary because of the geometrical and material nonlinearities of 
the interface material. The objective of this section is to obtain the tangent stiffness matrix K® and the 
internal force vector fj® t required in the nonlinear solution procedure. 

A 2n-noded isoparametric interface element with 6n degrees of freedom and applicable to three dimen- 
sional analysis is used. The element consists of an upper and lower surface Sf 1 with n-nodes each. The 
natural coordinate system is rji and 772 - For the surfaces S^, node j has three translational degrees of free- 
dom qfj , q^j , q%j with the first subscript implying the associated global direction. The nodal displacement 
vector q is arranged as follows, 

(4.1) q = {q + ,q~} T 

^ = tifv 9si. 9 *.. vinY r 

The displacement field 772 ), j = 1,2,3 of the upper and lower surfaces of the element is arranged in 

the vector UA The displacement vector U ;l is related to the global displacement degrees of freedom vector 
q 1 * 1 as, 

(4.2) U ± =Nq ± 


where N is a 3 X 3n shape function matrix. Substituting Equation 4.2 into Equation 3.11 gives the virtual 
relative displacement vector 67\ in terms of (^q^ 

(4.3) = [R t N (5q+ , -R t N <5q“] 

It is convenient to define the 3x3 n matrix B s , 

(4.4) B s = R t N 


and the 3x6 n relative-displacement/displacement matrix B, 

(4.5) 6%= [B s ,— B s ](5q = B(5q 


The internal force vector of the interface element is obtained by substituting Equation 4.5 in 3.1, 

(4.6) 5Wf nt = 4q r J j. B t T d,S™ = bq T f/;, 

where T is the traction vector acting on the deformed mid-surface and the integration is performed over 
the mid-surface of the deformed element. The tangent stiffness matrix stems from the linearization of the 
internal force vector and is obtained using Taylor’s series expansion about the approximation qW 


(4.7) 


fj e nt (q«+Aq)«C t (q«) + 


<9q 


Aq- 


J q=q(*) 
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where i is the iteration number. 

In numerical analyses, the internal force vector needs to be computed accurately, and the tangent stiff- 
ness matrix may be computed approximately. The computation of the tangent stiffness matrix is intensive 
and a very accurate expression is not required. Therefore, the derivative of the relative-displacement ma- 
trix/displacement matrix with respect to the nodal displacements are neglected. In addition, the partial 
derivatives of M in Equation 3.12 with respect to q are neglected. Thus, from Equation 4.7, the approxi- 
mate 6n x 6n tangent stiffness matrix is 


(4.8) 


B s = B+ = B 
B = [B S ,-B S ] 


where B implies approximation to B. 

(4.9) K e t = rj J J B t DB dS™ 

where D is the material tangent stiffness, and is later defined in detail. Equation 4.9 is rewritten using the 
relation in Equation 4.8, 


(4.10) 


Kf 


K s — K s 
-K s K s 


where K g is a 3n X 3 n submatrix defined as 


(4.11) 


K, = 


II, 


B s t DB b 


ds: 


The approximations for the tangent stiffness matrix save computational time. 


4.1. Material Tangent Stiffness. The components of the material tangent stiffness D are obtained 
in the incremental form, 

f)T 

(4.12) d'l) = = DijSAj 

First consider the case in which there is no interpenetration, that is, for A3 > 0. The components of D are 
obtained by differentiation of Equation 3.25 according to Equation 4.12, 


(4.13) 


D, 


rpc 

± i 

A c . 


, A i A J 

Sii - 


W /i' 


OL — (3 


I A? 


I ok — 2 


Q 


where 5{j is the Kronecker delta, Q is given by Equation 3.22, and w is defined by, 


(4.14) 


w = { 1 if d = fjPd, if d> fiP 


Now consider the case for which interpenetration is detected, that is, A3 < 0. The non-zero components 
of D are given by Equation 4.13 for i,j = 1, 2 and the component related to interpenetration, 

(4.15) H33 = K 0 ^1 + k |A 3 |^ exp 

where K 0 = T£ exp(l//?)/A§. The range of the values of D 33 should be restricted by two conditions: (1) A 
small D 33 induces interpenetration, and (2) a large D 33 produces ill-conditioned matrices. A list of references 
on these restrictions is given by Davila et al[13] . The value of D 33 should be in the range, 

(4.16) 10 6 N/mm 3 < D 33 < 10' N/mm 3 
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The upper bound of the condition cannot be guaranteed because of the exponential nature of Equation 4.15. 
Therefore, for A 3 < 0, the expressions T 3 and ID 33 are modified to have the form 

(4.17) T 3 = K 0 A 3> £>33 = K 0 

and Ko = T§ exp(l//3)/A§. 

The material tangent stiffness is non-symmetric, and can be positive definite, semi-definite, or negative 
definite. For fi > 1, the matrix Dij is negative definite. The material tangent stiffness matrix has properties 
of an anisotropic material, one which has strong dependence on the relative displacements in all directions. 
For single-mode delamination, D is fully diagonal, otherwise, some of the off-diagonals are non-zero. 

4 . 2 . Consistent and Inconsistent Tangent Stiffness. For the full-Newton-Raphson nonlinear so- 
lution procedure, the consistent tangent stiffness matrix is used in the finite element analysis. However, 
when softening constitutive laws with the consistent tangent stiffness are employed, the tangent stiffness 
matrix is often ill-conditioned and a converged solution may not be obtained[22]. An alternative solution is 
to refine the mesh ahead of the crack-tip or to decrease the maximum interfacial strength[15, 23]. Refining 
the mesh size increases the computational time, and lowering the maximum interface strength can result in 
a premature initiation of delamination [7]. Alternatively, researchers often utilize a positive definite matrix 
such as the material secant stiffness when dealing with softening constitutive laws. However, a large number 
of iterations results in using the material secant stiffness. As an alternative, three different modifications 
to the tangent stiffness matrix eliminate these convergence difficulties while a converged solution can be 
obtained in a small number of iterations: 

1. Equation 4.8, Kf d = max(0, Kf { ), i = 1,2, ..., 2 n 

2. Equation 4.9, K Sii = max(0, K Sii ), i = 1,2, ...,n 

3. Equation 4.13, Du = max(0, Du), i = 1, 2, 3 

The convergence rate of option 1 is better than option 2, and the convergence rate of option 2 is better than 
option 3. If the mesh is coarse, it is better to choose option 3. 

4 . 3 . Contact Elements. Interface elements were developed to model initial delaminated surfaces. All 
the components of the material tangent stiffness is zero, except for the case in which interpenetration is 
detected. If interpenetration is detected Equation 4.17 is used. Thus, these interface elements act like 
contact elements. 

5. Finite Element Results. Numerical results are presented for quasi-static loading and unloading of 
the double cantilever beam (DCB), the end load split (ELS), end notch flexure (ENF), and fixed ratio mixed 
mode (FRMM) fracture test specimens. Results are also presented for quasi-static loading of the mixed 
mode bending (MMB). Mode I fracture occurs in the DCB, Mode II occurs in the ELS and ENF, and Mode 
I and II occur in the FRMM and MMB. The fracture test specimen configurations are shown in Figure 6.5. 

Mode I and mixed-mode test specimens are modeled with the laminate stacking sequence [0]]] and 
the unidirectional material properties of Graphite- Epoxy listed in Table 6.1. An isotropic material with 
E = En and v = V 12 are used for the Mode II test specimens rather than composite. The maximum 
interfacial strength and the critical energy release rates are listed in Table 6.2. The geometrical properties 
are the length L = 100 mm, the arm thickness h = 1.5 mm, and width B = 10 mm. For the DCB, the 
geometrical properties are different from the other test specimen configurations: L = 150 mm, h = 1.5 mm, 
and B = 20 mm. The initial crack length ao of each test specimen is: DCB - 50 mm, ENF - 30 mm, ELS - 
50 mm, FRMM - 40 mm, and MMB - 20 mm. 
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The interface elements are positioned between the upper 0° lamina and the lower 0° lamina. Delami- 
nation is constrained to grow in the plane between the upper and lower laminates. Interface elements with 
contact properties were placed along the initial crack length and interface elements formulated with the 
softening law are placed along the bonded length. The upper and lower laminates are modeled with C3D8I 
incompatible- mode 8 node solid element available in ABAQUS. Each lamina is modeled with one element 
through the thickness, 100 elements along the length of the laminate, and one element across the width. See 
Figure 6.6a. For the DCB, three elements along the width are used. The eight node isoparametric interface 
element for three-dimensional analysis shown in Figure 6.6b is compatible with C3D8I solid element. The 
element was implemented in the commercial finite element code ABAQUS as an UEL subroutine. Three 
point Gauss integration is used for the computation of the tangent stiffness matrix and internal force vector. 

An incremental-iterative approach is adopted for the nonlinear finite element analysis, and the Newton’s 
method available in ABAQUS is used to trace the loading path of the specimens with a displacement-control 
analysis. For the MMB, the Riks method available in ABAQUS is used. The modification to the tangent 
stiffness matrix used is option 2 discussed in the section of interface elements. The response of the test 
specimens is characterized by the load-deflection curve. A typical finite element model of one of the test 
specimens consists of about 300 elements, and 2000 degrees of freedom. The computational time required 
was about 1200 seconds of CPU time on a Sun Solaris 2000. The average number of iterations per load 
increment is 7. 

The finite element solutions are compared to the beam analytical solutions derived from linear elastic 
fracture mechanics. The analytical solutions for the DCB and ENF are given by Mi et al. [15] , and for the 
FRMM and ELS are given by Chen et al. [23] . The finite element solutions for the MMB test specimen are 
compared to the analytical solution in the appendix. 

The DCB shown in Figure 6.5a is used to determine the interlaminar fracture toughness in Mode I. The 
displacement w is symmetrically applied, equal and opposite at the tip of the upper and lower arm of the DCB 
test specimen. The corresponding reaction force P is computed. The other end of the specimen is clamped. 
The response of the DCB is shown in Figure 6.7a. For a loading-unloading cycle, excellent agreement of the 
FEM results are obtained compared to to the closed form solutions and to the experimental data([13]). The 
contour plot of the effective von Mises stress of the DCB is shown from the top view near the delamination 
front in Figure 6.7b. The black strip is a region of low stress values, indicating that delamination has 
occurred. The gray strip is a region of intermediate stress values, and is where the material points are 
softening. The boundary of the black and gray strip is the location of the crack-tip. The white strip is the 
location of high stresses, and is the region where onset of delamination is occurring. Non-self-similar crack 
growth occurs because of the anticlastic bending effect. The tangent stiffness matrix in the Newton-Raphson 
methods did not converge at the limit point because of the large value of the maximum interfacial strength, 
T c . Any of the modifications to the tangent stiffness matrix discussed in the section of interface elements, 
produced converged solutions. 

The ELS and ENF are shown in Figure 6.5b and 6.5c are used to determine the interlaminar fracture 
toughness in Mode II. For the ELS, the load P is applied at the tip such that the lower arm of the ELS 
remains in contact with upper arm. The other end of the specimen is clamped. The ENF specimen is simply 
supported, and the downward vertical displacement is specified at the mid-span of the specimen. The 
corresponding reaction force P 2 is computed. The response of the ELS and ENF are shown in Figure 6.8a 
and 6.8b, respectively. For a loading-unloading cycle, excellent agreement of the FEM results are obtained 
compared to the closed form solutions. 
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The FRMM is shown in Figure 6.5d, and is used to evaluate empirical failure criteria for mixed-mode 
delamination. The displacement w is specified at the tip of the upper arm and the corresponding reaction 
force P is computed. Mode I is 43% and Mode II is 57%. The response for a = 2 and a = 4 (see Equation 
2.1 and 2.2) is shown in Figure 6.9a and 6.9b respectively. For a loading-unloading cycle, excellent agreement 
of the FEM results are obtained compared to the closed form solutions. 

The MMB is shown in Figure 6.5c, and is used to evaluate empirical failure criteria for mixed-mode 
delamination. The length of the lever arm c described in the report by Reeder [20] is chosen such that 
the mixed mode ratio from pure Mode I to pure Mode II can be varied. In this paper, c = 43.72 mm, 
so that the Mode I and Mode II contributions are 50% each. The MMB is simply supported, and two 
proportional loads are applied. The load Pi is applied upward at the tip of the upper arm, and another load 
P 2 is applied downward at the mid-span. During the loading, the ratio P 1 /P 2 = 2c/(2c + L) is fixed. The 
responses for a = 4 and a = 2 are shown in Figure 6.10a and 6.10b. The finite element response is compared 
to the analytical solutions in the appendix. In the first analysis, geometric nonlinearity is used. In the 
second analysis both geometric linearity and nonlinearity are compared with the analytical solutions. The 
discrepancies on the response corresponding to the stable crack growth of the load-deflection response are 
because the analytical solution does not consider the effects of geometric nonlinearities. Excellent agreement 
is obtained with the analytical solutions. 

6. Concluding Remarks. The mechanics of an idealized interface material is presented. The kinemat- 
ical relations and the irreversible constitutive law mathematically describe the mechanics of the delamination 
process. The delamination process comprises three interrelated phases: the initiation of delamination, the 
evolution of the degradation zone, and the delamination growth. It predicts initiation of delamination based 
on a multi-axial stress criterion, and progression of delamination based on an empirical mixed-mode fracture 
criterion. A damage parameter is included to prevent the restoration of the previous cohesive state between 
the interfacial surfaces. The constitutive law is implemented with interface elements using the principal of 
virtual work. To demonstrate the irreversibility capability of the constitutive law, steady-state crack growth 
is simulated for quasi-static loading-unloading cycle of various fracture test specimen configurations. The 
finite element solutions are in excellent agreement with the analytical solutions. 
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Fig. 6.5. Fracture test specimen configurations. 
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Fig. 6.6. Finite element modeling. 
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(a) Load-opening response of the DCB (b) Non-self similar delamination growth 


Fig. 6.7. DCB with ao — 50mm 
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FlG. 6.8. Load- displacement response of Mode II test specimen configurations. 
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Fig. 6.9. Load- displacement response of the FRMM test specimen configurations. 
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Fig. 6.10. Load-displacement response of the MMB test specimen configurations. 
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Appendix A. Mixed-Mode Analytical Solutions. The beam analytical solutions based on linear 
elastic fracture mechanics for the MMB test specimen are presented without details. In general, the total 
energy release rate is 

(A.l) Gy = G[ + Gjj 

G / and G / / are the Mode I and Mode II energy release rate contributions. The delamination propagates 
when, 


(A.2) 


/-~r s~im . /~im 

Lit — W ^ I I ^ II 


and G c is the critical energy release rate, G™ and G™ are the the Mode I and Mode II energy release rates 
at crack propagation. For all the fracture test specimens, it is possible to express <j> = G™/G/}, where 
4> G [0, oo), so that the G c value can be computed based on the fracture criterion in Equation 2.2 


(A.3) 


G c 


(1 + 0 ) 




-( 2 /«) 


The derivations to obtain the expression of <fi for the MMB specimen are omitted here, and is 


(A- 4) 


G/ Gf 4 / 6c — A \ 
G/7 ~~ G/} _ 3 \2c + Lj 


where c is the length of the lever arm in the test fixture ([20]), and A is the length of the MMB specimen. 
For simplifying purposes, the loads P/ and P/j associated to modes 1 and II respectively are defined as 


(A.5) 


Pi = - 

c 


6c -L 
4 A 


Pn = - 

c 


2 c+ L 


The load Pi is defined in Figure 6.5c. The initial load-deflection response is linear and given by 


(A.6) 


w i = 


16i /6c -L 


3c 


4 L 


Piaj 

El 


where E is the Young’s Modulus and I is the moment of inertia. The load-deflection response, when 
delamination propagates with a < i/2 is 


(A.7) 


16 Pi f 8BEIG C \ 3/2 
3 El \64P/ + 3 Pfj ) 


where B is the width of the beam. The load-displecement relation when delamination propagates with 
a > i/2 is obtained by solving the quadratic equation for a, 


(A. 8) (64 Pj + 3 P?! - 64 P/P//) a 2 

- (6PJ/P - 32P/P//A) a + (SP^i 2 - 8 BEIG C ) = 0 


and substituting its solution into 
(A.9) 


wi 


16 i /6c — A \ Pi a 3 
~37 / 4 A ) El 
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